Method For Predicting Lithology And Porosity From Seismic Reflection Data

ABSTRACT

Method for predicting lithology and porosity of subsurface rocks from seismic reflection data. The seismic data is inverted to yield elastic properties of the rocks such as the compressional and shear impedances. A rock physics model is built to relate porosity, the shale volume fraction, the fluid content of the rock and the elastic properties of the rock. The model is run backward in a second inversion process to solve for porosity and lithologic properties such as the shale volume fraction.

This application claims the benefit of U.S. Provisional Application No.60/574,901, filed May 27, 2004.

FIELD OF THE INVENTION

This invention relates generally to the field of characterizinghydrocarbon reservoirs and, more particularly, to lithology inversionand methods for predicting sand (or shale) volume and porosity insub-surface rocks. Specifically, the invention is a method forpredicting lithology and porosity from seismic reflection data.

BACKGROUND OF THE INVENTION

In characterizing hydrocarbon reservoirs, estimating reserves, anddeveloping models for how to best extract the hydrocarbons, it is usefulto know the lithology (for example, relative amounts of shale and sand)and associated porosity of the rocks in the target interval. Rockproperties can be measured directly from rock samples obtained fromwells but such samples are generally very limited in availability due tothe expense of drilling those wells. These properties can also beinferred from seismic data. Because of the complicated nature of thetheoretical relationships between the seismic data (reflectivity) andthe important rock properties (lithology, porosity, and fluid content),these two quantities are often related in practice through empiricalrelationships derived at wells, where both seismic and well measurementscoexist. These empirical relationships are then applied to the entirevolume of seismic data (or attributes derived from them) in order tomake predictions about rock properties away from the wells. The problemis that empirical models require a statistically significant sampling ofdata and yet the wells provide very limited and generally biased samplesof the reservoir properties. In regions where a large number of wellshave been drilled, pattern-based recognition methods and simpleempirical relationships can be used successfully to infer rockproperties from seismic data. However, in regions of limited wellcontrol, it is difficult to make accurate lithology predictions usingempirical relationships derived from just a few wells.

A commonly used method for determining clay content and porosity fromseismic data (or attributes of the seismic data) is to use linearregression to solve an equation of the following form:Impedance=A·φ+B·vshale+Cwhere φ is the porosity, vshale is the shale volume and A, B and C arethe constants that relate the porosity, vshale and impedances (or someother seismic attribute of interest) to one another. Regression methodsare more robust when they are used with larger datasets obtained fromwells penetrating different sections of the reservoir so that there is astatistically significant sampling of the data. In regions of limitedwell control the relationships derived in this manner cannot be usedwith confidence.

Another class of methods used to predict clay content and porosity fromseismic data uses pattern recognition, often implemented with neuralnetworks, to construct the necessary relationships. These methods use atraining set to identify patterns between the well and the seismic dataand then classify the remainder of the seismic data set according to thepatterns observed in the training set. The resulting relationships canbe quite complicated (and certainly allow more complexity than thesimple linear regression of the above-stated equation), but they arestill fundamentally empirical relationships based on observations at thewell rather than on a physical description. Consequently, these methodssuffer from the same problem as the regression methods in that theyrequire enough data examples (wells) in order to train the networkcompetently. With sufficient well control they can be very goodinterpolators (although generally poor extrapolators). In regions oflimited of well control, they are unreliable interpolators (as well aspoor extrapolators).

SUMMARY OF THE INVENTION

In one embodiment, the invention is a method for predicting lithologicproperties and porosity of a subsurface formation from seismic data,comprising:

-   -   (a) inverting the seismic data to obtain one or more bulk        elastic properties of the subsurface formation;    -   (b) constructing a rock physics model of the subterranean        formation, said model relating the lithologic properties,        porosity and fluid content to the bulk elastic properties of the        formation rock, said model comprising the following two        features: (i) compliances and densities of sand and clay mineral        fractions of the rock are characterized independently with        separate pore spaces, different pore aspect ratios, and        potentially different fluid types, and (ii) effective bulk and        shear elastic moduli are computed using a combination of        differential effective medium theory and Gassman fluid        substitution;    -   (c) building a fluid fill model indicating the type of fluid        present at each location in the subsurface;    -   (d) computing in tabular form values of said one or more elastic        properties as predicted by the rock physics model for a range of        possible values for said porosity and lithology properties in        each fluid type present in the model and then numerically        computing corresponding tables of the derivatives of the elastic        properties with respect to porosity and clay content; and    -   (e) using the computed tables of the elastic properties and        their derivatives, along with the fluid type information to        minimize a pre-selected objective function and thereby invert        the rock physics model to obtain the lithologic properties and        the porosity from the bulk elastic properties and fluid content        information for the formation.

Typical bulk elastic properties include compressional impedance, shearimpedance, bulk modulus, shear modulus, compressional velocity, shearvelocity or any other elastic parameters. Typical lithologic propertiesinclude the volume fractions of shale (clay) and of sand.

In some embodiments of the invention, the rock physics model has a solidmatrix composed of sands and clays and a total pore space partitionedinto clay-related pores and sand-related pores, and the clay-relatedpores are assumed to be filled primarily by water (actually brine). Insome of these embodiments, the pressure equalization assumptions for theclay pores differ from the sand pores. For example, the brine-filledclay-related pores may be added during the effective media calculationso that only the sand-related pores are filled using Gassman theory.This corresponds to the mixed frequency case where the pressure withinthe clay pores is not equalized during the passage of a seismic wavewhile the pressure in the larger sand pores is equalized. In otherembodiments of the present invention, both the sand-related andclay-related pores are empty during the differential effective mediumcomputation and are later filled with fluid using Gassman theory.

In some embodiments of the invention, the second inverting step solvesfor the lithologic properties and porosity using an iterative processand converging to a solution by minimizing the squared difference, oroptimizing the L₁ norm of the difference, between the bulk elasticproperties obtained from the seismic data and the values obtained forthe same properties by forward modeling with the rock-physics model. Inother related embodiments, the iterative process converges to a solutionby finding a maximum a posteriori estimate (MAP) of the lithologicproperties and porosity using model and data covariance matricesestimated from well data and inversion results at such well. In some ofthe above-described embodiments, the iterative process is Newton-Raphsoniteration. To speed up and simplify the solution process, someembodiments of the invention further comprise constructing a table of Pand S impedances that have been forward modeled with the rock physicsmodel for representative values of clay content and porosity; andpre-computing tables of the derivatives of the P and S impedances withrespect to porosity and clay content using finite differenceapproximations. The second inverting step is then performed by takingthe pre-computed tables and performing a non-linear inversion todetermine the combination of clay content and porosity that isconsistent with the P and S impedances at each point in the seismic datavolume.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention and its advantages will be better understood byreferring to the following detailed description and the attacheddrawings in which:

FIG. 1 is a flow chart showing the primary steps of one embodiment ofthe present invention; and

FIG. 2 is a three-dimensional view of sand bodies as predicted by thepresent inventive method from actual seismic data.

The invention will be described in connection with its preferredembodiments. However, to the extent that the following detaileddescription is specific to a particular embodiment or a particular useof the invention, this is intended to be illustrative only, and is notto be construed as limiting the scope of the invention. On the contrary,it is intended to cover all alternatives, modifications and equivalentsthat may be included within the spirit and scope of the invention, asdefined by the appended claims.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The present invention is a method for inferring the clay content andporosity of a reservoir from seismic reflection data. [The terms “shale”and “clay” are used interchangeably herein.] It is based upon thephysics of wave propagation through elastic media and a rock physicsmodel relating the elastic properties of rocks to their grain and fluidcomponents and their micro-pore structure. Therefore it mathematicallyrelates the lithologic description of the rocks to their seismicresponse and does not rely on empirical models. Because a theoreticalmodel of the rock physics is used, a statistically significant samplingof “ground truth” is not needed, and the method can be applied withlimited well control early in the development cycle when uncertaintiesare highest and economic risk is greatest. The rock porosity andlithology are predicted simultaneously, ensuring that these twoquantities are consistent with both the data and the rock physics model.Furthermore, this approach can simultaneously predict the reservoirquality in both the hydrocarbon and brine legs of a reservoir and doesnot require separate calibrations in the various fluid phases as isnecessary with empirical methods.

The method consists of three steps, as illustrated in FIG. 1. In step 1,the seismic data are inverted for elastic properties such ascompressional (P) and shear (S) impedances (I_(p) and I_(s),respectively) using standard techniques. (The seismic data are shownseparated into at least two constant angle stacks 4 because this isnecessary in order to reliably obtain both the P and S impedances.)Other elastic properties such as the bulk modulus and shear modulus, orcombinations thereof such as compressional and shear velocity may alsobe used. Such inversion techniques are described in, for example, T.Tonellot, D. Mace, V. Richard, “Prestack elastic waveform inversionusing a priori information,” SEG Expanded Abstracts (1999); James J.Carazzone and Leonard J. Srnka, “Elastic Inversion of Gulf of MexicoData, in Offset-Dependent Reflectivity”, Theory and Practice of AVOAnalysis, edited by John P. Castagna and Milo M. Backus, SEG (1993);Arild Buland, Martin Landro, Mona Andersen, and Terje Dahl, “AVOInversion of Troll Field data,” Geophysics, 1589-1602 (1996). In Step 2,the fluid type (gas, oil, brine, etc.) is defined for each and everypoint in the elastic properties volume and in this manner a fluid fillmodel is built. In step 3, the elastic properties obtained in step 1 arecombined with fluid information compiled in step 2 and converted tolithology and porosity (fractional shale volume V_(sh) and porosity φare shown in FIG. 1) using an appropriate rock physics model. Step 3 isapplied as a second or cascaded inversion following the seismicinversion in step 1.

In order to perform the second inversion, the present invention uses arock physics model that relates porosity, volume of shale or clay,V_(SH), and fluid content to the bulk elastic properties of the rocksuch as P-impedance and S-impedance. A preferred model in clasticenvironments is the shaly-sand mixture model described in Xu and White,“A new velocity model for clay-sand mixtures,” Geophysical Prospecting43, 91-118 (1995), and in Xu and White, “A physical model for shear wavevelocity prediction,” Geophysical Prospecting 44, 687-717 (1996), ormodifications of that model such as those described herein. The Xu-Whitemodel is complex, and inverting it poses a major problem. The presentinventive method solves this problem in a practical and efficient way.

The model has two key features. First, the compliances of the sand andclay mineral fractions of the rock are characterized independently withseparate pore spaces and different effective pore aspect ratios. Second,the bulk and shear elastic moduli of the dry frame are computed using acombination of the scattering theory of Kuster and Toksöz's “Velocityand attenuation of seismic waves in two-phase media: Part 1: Theoreticalformulation”, Geophysics 39, 587-606 (1974) and the differentialeffective medium theories of Bruner in “Comment on Seismic Velocities inDry and Saturated Cracked Solids by Richard J., O'Connell and BernardBudianskey”, Journal of Geophysical Research 81, 2573-2576, (1976) andCheng and Toksöz in “Inversion of seismic velocities for pore aspectratio spectrum of a rock”, Journal of Geophysical Research, vol 84, pp.7533-7543 (1979). Equations disclosed by Gassman in “Elasticity ofporous media”, Vierteljahrschift der Naturforschenden in Zürich, vol 96,pp. 1-21 are then applied to obtain the low frequency velocity in thefluid-saturated rock. This model computes relationships betweenvelocity, density, clay content and porosity that are explicit,consistent and physically based. As a result, a large number of nearbywells or assumed analogs are not required to characterize thesubsurface. The model is next described in more detail for certainembodiments of the present invention.

The mathematical expressions in the selected rock physics model providea method for determining P and S velocities and densities in rocks givenclay content, porosity, water saturation and fluid properties. Theseparameters can be recombined to give the impedances or any other set ofisotropic elastic properties that are produced in step 1. A typicalmodel assumes a solid matrix made of sands and clays. The total porespace can be partitioned into clay-related pores and sand related pores.If φ denotes total porosity, thenφ=φ_(S)+φ_(C)  (1)where φ_(S) is the portion of the rock occupied by stiff or sandstonepores, and φ_(C) is the porosity associated with compliant or shalepores. Fractional shale volume V_(SH) and fractional sand volume V_(SS)are used to estimate φ_(C) and φ_(S). Since V_(SH)+V_(SS)+φ=1, assumingthat φ_(C) and φ_(S) are proportional to V_(SH) and V_(SS),respectively, implies that $\begin{matrix}{\phi_{C} = {V_{SH}\frac{\phi}{1 - \phi}}} & (2) \\{and} & \quad \\{\phi_{S} = {V_{SS}{\frac{\phi}{1 - \phi}.}}} & (3)\end{matrix}$Having divided the pore space into compliant and stiff pores, the effectof pore shape on the elastic properties of the composite can beestimated using the following equations from the 1974 Kuster and Toksözpaper: $\begin{matrix}{{K_{0} - K_{m}} = {\frac{1}{3}( {K^{\prime} - K_{m}} )\frac{{3\quad K_{0}} + {4\quad\mu_{m}}}{{3\quad K_{m}} + {4\quad\mu_{m}}}\phi_{C}{T_{iijj}( \alpha_{C} )}}} & (4) \\{and} & \quad \\{{\mu_{0} - \mu_{m}} = {\frac{( {\mu^{\prime} - \mu_{m}} )}{5}\frac{\begin{matrix}{{6\quad{\mu_{0}( {K_{m} + {2\quad\mu_{m}}} )}} +} \\{\mu_{m}( {{9\quad K_{m}} + {8\quad\mu_{m}}} )}\end{matrix}}{5\quad{\mu_{m}( {{3\quad K_{m}} + {4\quad\mu_{m}}} )}}\phi_{C}{F( \alpha_{C} )}}} & (5) \\{where} & \quad \\{{F(\alpha)} = {{T_{ijij}(\alpha)} - \frac{T_{iijj}(\alpha)}{3}}} & (6)\end{matrix}$in which K₀, K_(m), and K′ are the bulk moduli of the effective mediumwith clay pores only, the rock matrix, and the pore inclusion material,respectively, and μ₀, μ_(m), and μ′ are the corresponding shear moduli.μ′ is always zero for any pore fluids. α_(C) is the aspect ratios forcompliant pores (clay pores); and T_(iijj) (α) and F(α) are pore aspectratio functions derived from the tensor T_(ijkl) that relates theuniform strain field at infinity to the strain field within an elasticellipsoidal inclusion. The moduli of the rock matrix are a mixture ofthe sand and clay grain moduli, mixed using the Voight-Ruess-Hillaverage according to their relative proportions as given by the V_(SH)value. $\begin{matrix}{K_{m} = {\frac{1}{2}\{ {{V_{sh}K_{sh}} + {( {1 - V_{sh}} )K_{ss}} + \frac{1}{\lbrack {\frac{V_{sh}}{K_{sh}} + \frac{( {1 - V_{sh}} )}{K_{ss}}} \rbrack}} \}}} & (7) \\{and} & \quad \\{\mu_{m} = {\frac{1}{2}\{ {{V_{sh}\mu_{sh}} + {( {1 - V_{sh}} )\mu_{ss}} + \frac{1}{\lbrack {\frac{V_{sh}}{\mu_{sh}} + \frac{( {1 - V_{sh}} )}{\mu_{ss}}} \rbrack}} \}}} & (8)\end{matrix}$The symbols K_(sh) and μ_(sh) in equations (7) and (8) are the bulk andshear moduli of the clay minerals. The symbols K_(ss) and μ_(ss) inequations (7) and (8) are the bulk and shear moduli of the sandminerals. The equations for the pore aspect ratio functions T_(iijj)(α)and F(α) in equation (6) are given in the Appendix of the 1995 Xu andWhite article.

In one embodiment of the invention, a key modification relative to thepublished models is the further assumption that clay pores are filledprimarily with bound water, because of the relatively high capillarypressure in clay pores as a result of the extremely small sizes of clayparticles. A further modification is to assume that the pressure inthese small pores is not equalized with the passage of a seismic wavebecause of their small pore throats. Mathematically, this is equivalentto saying they are unrelaxed and comprise a high-frequency component ofthe system. Since sand-related pores tend to be much bigger thanclay-related pores, pore fluids in sand-related pores can be equalizedeasily. At seismic frequencies, these pressure-equalized sand-relatedpores can be treated mathematically as being relaxed and comprising alow frequency component of the system. In this embodiment of the presentinvention, the Kuster-Toksöz equations are again used to calculate theelastic properties of the “dry rock frame” (sand pores only) letting K′and μ′ be zero. $\begin{matrix}{{K_{d} - K_{0}} = {{- \frac{1}{3}}K_{0}\frac{{3\quad K_{d}} + {4\quad\mu_{0}}}{{3\quad K_{0}} + {4\quad\mu_{0}}}\phi_{S}{T_{iijj}( \alpha_{S} )}}} & (9) \\{and} & \quad \\{{\mu_{d} - \mu_{0}} = {{- \frac{\mu_{0}}{5}}\quad\frac{{6\quad\mu_{d}( {K_{0} + {2\quad\mu_{0}}} )} + {\mu_{0}( {{9\quad K_{0}} + {8\quad\mu_{0}}} )}}{5\quad{\mu_{0}( {{3\quad K_{0}} + {4\quad\mu_{0}}} )}}\phi_{S}{{F( \alpha_{S} )}.}}} & (10)\end{matrix}$Here K_(d) and μ_(d) are the bulk and shear moduli of the “dry rockframe”, and α_(s) is the aspect ratios for stiff pores (sand pores).

However, the Kuster-Toksöz equations require φ/α<<1. Typical values foraspect ratios are 0.035 for shale pores and 0.12 for sandstone pores.Therefore, the Kuster-Toksöz equations are applicable only for very lowporosity. The differential effective medium method may be incorporatedinto the Kuster-Toksöz formulations to overcome this restriction. Toapply the differential effective medium method, the total porosity ispreferably modified using the following equation before partitioning thepore space.φ′=−ln(1−φ)  (11)The modified total pore space is then partitioned into sets of pores sothat the pore concentration for each set satisfies the Kuster-Toksözcondition. Beginning with solid rock, the Kuster-Toksöz equations areused to compute the effective medium that results from adding a smallset of pores to the matrix. In another key modification relative to themethod described in Xu and White (1995), the small set of pore space isthen divided into sand-related and clay-related portions using equations(1) to (3). The Kuster-Toksöz equations (4) to (10) are then used tocalculate the effect of clay-related and sand-related pores on elasticproperties. The process is repeated, using the effective medium from theprevious calculation as the rock matrix for the next calculation, untilthe total pore volume has been added to the rock. Finally, Gassmann'sequation (referenced previously) is used to put pore fluids into thesand pores. $\begin{matrix}{{K = {K_{d} + \frac{( {1 - \frac{K_{d}}{K_{0}}} )^{2}}{\frac{\phi_{S}}{K_{f}} + \frac{( {1 - \phi_{S}} )}{K_{0}} - \frac{K_{d}}{K_{0}^{2}}}}},} & (12)\end{matrix}$  μ=μ_(d),  (13)ρ=φρ_(f)+(1−φ)ρ₀.  (14)After obtaining the effective bulk and shear moduli, P- and S-wavevelocities can be calculated using the following equations:$\begin{matrix}{{V_{P} = \sqrt{\frac{( {K + {\frac{4}{3}\mu}} )}{\rho}}},} & (15) \\{and} & \quad \\{V_{S} = {\sqrt{\frac{\mu}{\rho}}.}} & (16)\end{matrix}$In step 3 of FIG. 1, the impedances found in the step 1 are inverted forporosity and V_(SH) for every point in the seismic volume using (in thedescribed embodiment) the Xu-White forward model or its modificationsdescribed above. The porosity and V_(SH) values that best fit theimpedances are found by minimizing an objective function. Typically, theobjective function will consist of a term measuring the match betweenthe bulk elastic properties generated by the rock physics model and thebulk elastic properties obtained from seismic inversion, and a termconstraining the predicted lithology and porosity. The simplestobjective function is the least-squares objective function containingthe squared difference between the forward modeled impedances and theobserved impedances.ε²=(I _(p) −XW _(p)(φ,V _(sh)))²+(I _(s) −XW _(s)(φ,V _(sh)))²  (17)Here, I_(p) and I_(s) are the impedances derived in the first step, andXW_(p) and XW_(s) are the forward modeled p and s impedances using theseries of equations (1-15). As previously mentioned, other sets ofisotropic elastic parameters could be used, depending on what wasproduced in step 1. The minimization of the squared error is anon-linear problem. A preferred method for solving this problem is aNewton-Raphson iteration (see W. Press, et al., Numerical Recipes: TheArt of Scientific Computing, Cambridge University Press (1986) pp.254-259). An initial guess at the solution is made, and then it isiteratively updated it by solving the linearized equation set$\begin{matrix}\begin{matrix}{{\begin{bmatrix}\frac{\partial{{XW}_{p}( {\phi_{k},V_{{sh}_{k}}} )}}{\partial\phi} & \frac{\partial{{XW}_{p}( {\phi_{k},V_{{sh}_{k}}} )}}{\partial V_{sh}} \\\frac{\partial{{XW}_{s}( {\phi_{k},V_{{sh}_{k}}} )}}{\partial\phi} & \frac{\partial{{XW}_{s}( {\phi_{k},V_{{sh}_{k}}} )}}{\partial V_{sh}}\end{bmatrix}\begin{bmatrix}{\Delta\quad\phi_{k}} \\{\Delta\quad V_{{sh}_{k}}}\end{bmatrix}} =} \\{\begin{bmatrix}{I_{p} - {{XW}_{p}( {\phi_{k},V_{{sh}_{k}}} )}} \\{I_{s} - {{XW}_{s}( {\phi_{k},V_{{sh}_{k}}} )}}\end{bmatrix}.}\end{matrix} & (18)\end{matrix}$The subscript k in equation (18) is an iteration index. In order toevaluate the coefficients on either side of equation (18), one canevaluate the modified Xu-White model at the current guess φ_(k),V_(sh)_(k) as well as the derivatives of the model with respect to shalevolume and porosity. The equations are then solved for a model updateΔφ_(k),ΔV_(sh) _(k) . This update is added to the current model.Iteration continues until the solution converges.

Equation set (18) must be solved at a substantial number of points inthe seismic volume. Evaluation of the modified Xu-White model usingequations (1)-(17) is time-consuming. Iteratively evaluating themodified Xu-White model for all points in a seismic volume would be verycomputationally intensive. Furthermore, the derivation of analyticexpressions for the derivatives on the left-hand side of equation (18)is not tractable. To circumvent these difficulties, one can construct atable of P and S impedances (velocities scaled by density) that havebeen forward-modeled for representative combinations of clay content (0to 100%) and porosity (0 to 40% in siliciclastics). Tables of thederivatives of the P and S impedances with respect to porosity and claycontent are also pre-computed using finite difference approximations.The inverse rock physics modeling (running the model backwards) is doneby taking the pre-computed tables and performing a non-linear inversionto determine the combination of clay content and porosity that isconsistent with the P and S impedances derived at each point in theseismic volume.

Equations (1)-(17) describing the modified Xu-White rock physics modeldepend on the properties of the fluid filling the pore space through thebulk modulus K_(f) and density ρ_(f) of the fluid filling the sand poresin equations (12) and (14) respectively, as well as through the bulkmodulus of the fluid filling the clay pores K′ in equation (4). Asdescribed previously, the clay pores are filled with brine. The sandpores however, are filled with the appropriate reservoir fluid, eitherbrine or a combination of brine and hydrocarbons. The brine andhydrocarbons are mixed in proportion to the water saturation. Typically,a fixed effective water saturation is specified. In order to apply thesecond inversion, the pore fluid at every point in the seismic volumemust be specified (step 2 in FIG. 1). This involves defining the regionsof the reservoir corresponding to the hydrocarbon leg (either gas and/oroil) and the brine leg. The regions of the seismic volume correspondingto the different fluids can be constructed from the interpreted seismichorizons and polygons defining the reservoir and from the depths of thefluid contacts either penetrated in the wells, or inferred from pressuremeasurements taken in multiple wells, or from geophysical evidence suchas direct hydrocarbon indicators. Each region is identified in thevolume with a fluid identification flag. This fluid identification flagis associated with bulk moduli and densities appropriate for the fluidtype (gas, oil, or brine) and its environmental (pressure andtemperature) conditions as well as its relevant compositional parameters(gas gravity for gas, API and gas oil ratio for oil, and salinity forbrine). All these factors are determined from relevant measurements madein wells, and then translated into fluid bulk moduli and densities usingappropriate models such as those described by Michael Batzle and ZhijingWang, “Seismic properties of pore fluids,” Geophysics, 1396-1408 (1992).The Xu-White tables and its derivatives used in the application ofequation (18) during the rock physics inversion must be pre-computed foreach unique fluid identification flag in the fluid ID volume. Theappropriate set of tables is then applied at each point in the volumeduring the rock physics inversion.

Some of the constants used in the rock physics modeling requirecalibration. In particular, the pore aspect ratios of the sand and clay,as well as the grain properties of the clay are preferably adjusted tofit available well data prior to performing the rock, physics inversionof the seismically derived elastic properties. This calibration stepconsists of using shale volume and porosity values derived from welllogs along with the appropriate fluid properties to forward-model theelastic properties and compare them with measured density and sonic logsin the well. The rock physics parameters are adjusted until a reasonablematch between predicted and recorded well logs is obtained.

Because the elastic properties derived from the inversion in step 1 ofFIG. 1 can be noisy and therefore inconsistent with the modeledimpedances, strictly minimizing a least squares objective function as inequation (17) can produce noisy estimates of the rock properties. It isoften desirable to modify the objective function so as to add someconstraints in the rock physics inversion (step 3). A common strategy isto perform damped least squares by adding an additional term to theobjective function in which the model with the minimum squared length issought. A number of other strategies can be employed, as detailed invarious textbooks on geophysical inversion such as Menke, Geophysicaldata analysis: discrete inverse theory, Academic Press (1984). Apreferred approach is to use the maximum likelihood estimator detailedin Tarantola, Inverse problem theory: methods for data fitting and modelparameter estimation, Elsevier Science Publishers (1987). Appropriatedata and model covariance matrices are derived from data measured at thewell locations during the calibration step. Equation (17) measures theerror using an L₂ norm. Other norms, such as an L₁ norm, can be used tomeasure the error as well.

EXAMPLES

The present inventive method was applied to seismic data acquired over areservoir containing gas, oil and brine. FIG. 2 shows a 3-D image of theinferred sand channel winding through the inverted V_(SH) volume (theshaly parts have been made invisible) produced by the present inventivemethod. Two wells drilled in the reservoir confirm the location of thesands and the volume percentages. In addition to aiding in wellplacement, accurate sand/shale volume predictions such as thatrepresented by FIG. 2 can be used to estimate reserves, both of whichare important in the early development phases of a field when wellcontrol is limited and business risk is greatest.

The foregoing description is directed to particular embodiments of thepresent invention for the purpose of illustrating it. It will beapparent, however, to one skilled in the art, that many modificationsand variations to the embodiments described herein are possible. Forexample, persons skilled in the art will know of other modifications tothe differential effective media model. All such modifications andvariations are intended to be within the scope of the present invention,as defined by the appended claims.

1. A method for predicting lithologic properties and porosity of asubsurface formation from seismic data, comprising: (a) inverting theseismic data to obtain one or more bulk elastic properties of thesubsurface formation; (b) constructing a rock physics model of thesubterranean formation, said model relating the lithologic properties,porosity and fluid content to the bulk elastic properties of theformation rock, said model comprising the following two features: (i)compliances and densities of sand and clay mineral fractions of the rockare characterized independently with separate pore spaces, differentpore aspect ratios, and potentially different fluid types, and (ii)effective bulk and shear elastic moduli are computed using a combinationof differential effective medium theory and Gassman fluid substitution;(c) building a fluid fill model indicating the type of fluid present ateach location in the subsurface; (d) computing in tabular form values ofsaid one or more elastic properties as predicted by the rock physicsmodel for a range of possible values for said porosity and lithologyproperties in each fluid type present in the model and then numericallycomputing corresponding tables of the derivatives of the elasticproperties with respect to porosity and clay content; and (e) using thecomputed tables of the elastic properties and their derivatives, alongwith the fluid type information to minimize a pre-selected objectivefunction and thereby invert the rock physics model to obtain thelithologic properties and the porosity from the bulk elastic propertiesand fluid content information for the formation.
 2. The method of claim1, wherein the bulk elastic properties are selected from the groupconsisting of compressional impedance, shear impedance, bulk modulus,shear modulus, compressional velocity, shear velocity and other elasticparameters.
 3. The method of claim 1, wherein the lithologic propertiescomprise the volume fractions of shale (clay) and of sand.
 4. The methodof claim 1, wherein said model has a solid matrix composed of sands andclays and a total pore space partitioned into clay-related pores andsand-related pores, said clay-related pores being assumed to be filledprimarily by bound water during the differential effective mediacomputation and only the sand-related pores are filled using Gassmantheory.
 5. The method of claim 1, wherein said model has a solid matrixcomposed of sands and clays and a total pore space partitioned intoclay-related pores and sand-related pores, and both are empty during thedifferential effective media computation and later filled with fluidusing Gassman theory.
 6. The method of claim 1, wherein the inversion ofthe rock physics model solves for the lithologic properties and porosityusing an iterative process and converging to a solution by minimizingthe squared difference between the bulk elastic properties obtained fromthe seismic data and the values obtained for the same properties byforward modeling with the rock-physics model.
 7. The method of claim 6,wherein the iterative process is Newton-Raphson iteration.
 8. The methodof claim 1, wherein the objective function is the least-squaresobjective function containing the squared difference between elasticproperties generated by the rock physics model and the elasticproperties obtained from the seismic inversion step.
 9. The method ofclaim 1, wherein the objective function is a damped least-squaresobjective function.
 10. The method of claim 1, wherein the inversion ofthe rock physics model solves for the lithologic properties and porosityusing an iterative process and coverging to a solution by optimizing theL₁ norm of the difference between bulk elastic properties obtained fromthe seismic data and values obtained for the same properties by forwardmodeling with the rock-physics model.
 11. The method of claim 1, whereinthe inversion of the rock physics model solves for the lithologicproperties and porosity using an iterative process and converging to asolution by finding a maximum a posteriori estimate (MAP) of thelithologic properties and porosity using model and data covariancematrices estimated from well data and inversion results at the well.